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1.1 Introduction 

The area of excited-state dynamics is receiving increasing attention for a num- 
ber of reasons, including the importance of photochemical processes in basic 
energy sciences, improved theoretical methods and the associated theoretical 
understanding of photochemical processes, and the advent of femtosecond 
(and now attosecond) spectroscopy allowing access to more detailed experi- 
mental information about photochemical processes. Since photophysical and 
chemical processes are more complex than thermal (i.e., ground state) pro- 
cesses, simulations quickly become expensive and even unmanageable as the 
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Fig. 1.1. Schematic representation of potential energy surfaces for piiotopiiysical 
and photochemical processes: So, ground singlet state; Si, lowest excited singlet 
state; Ti, lowest triplet state; ABS, absorption; FLUO, fluorescence; PHOS, phos- 
phorescence; ISX, intersystem crossing; CX, conical intersection; TS, transition 
state. 



model system becomes increasingly realistic. With its combination of simplic- 
ity and yet relatively good accuracy, TDDFT has been finding an increas- 
ingly important role to play in this rapidly developing field. After reviewing 
some basic ideas from photophysics and photochemistry, this chapter will 
cover some of the strengths and weaknesses of TDDFT for modeling pho- 
toprocesses. The emphasis will be on going beyond the Born-Oppenheimer 
approximation. 

There are distinct differences between how solid-state physicists and chem- 
ical physicists view photoprocesscis. We believe that some of this is due to 
fundamental differences in the underlying phenomena being studied but that 
much is due to the use of different approximations and the associated lan- 
guage. Ultimately anyone who wants to work at the nanointerface between 
molecules and solids must come to terms with these differences, but that is 
not our objective here. Instead we will adapt the point of view of a chemical 
physicist (or physical chemist) see e.g., Ref. [Michl 1990]. 

The usual way to think about molecular dynamics is in terms of the po- 
tential energy surfaces that come out of the Born-Oppenheimer separation. 
In thermal processes, vibrations arc associated with small motions around 
potential energy surface minima. Chemical reactions are usually described 
as going over passes (transition states) on these hypersurfaces as the system 
moves from one valley (reactants) to another (products). Photoprocesses are 
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much more complicated fFig. [LIT) . Traditionally they include not only pro- 
cesss that begin by absorption of a photon, but also any process involving 
electronically excited states, such as chemiluminescence (e.g., in fireflies and 
glow worms) where the initial excitation energy is provided by a chemical 
reaction. The Frank-Condon approximation tells us that the initial absorp- 
tion of a photon will take us from one vibronic state to another without 
much change of molecular geometry, thus defining a Frank- Condon region on 
the excited-state potential energy surfaces. The molecule can return to the 
ground state by emitting a photon of the initial wave length or, depending 
upon vibronic coupling and perturbations from surrounding molecules, the 
molecule may undergo radiationless relaxation to a lower energy excited state 
before emitting or it may even decay all the way to the ground state without 
emitting. If emission takes place from a long-lived excited state of the same 
spin as the ground state, then we speak of fiuorescence. If emission takes 
place from an excited-state with a different spin due to intersystem crossing, 
then we speak of phosphorescence. If it is unsure whether the emission is 
fluorescence or phosphorescence, then we just say the molecule luminesces. 
Because of the large variety of de-excitation processes, excited molecules usu- 
ally return too quickly to their ground state for the molecular geometry to 
change much. We then speak of a photophysical process because no chemical 
reaction has taken place. Thus fluorescence is usually described as an excited 
molecule relaxing slightly to a nearby minimum on the excited-state potential 
energy surface where it is momentarily trapped before it emits to the ground 
state. It follows that the photon emitted during fluorescence is Stokes shifted 
to a lower energy than the photon initially absorbed. 

Photochemical reactions occur when the excited molecule decays to a 
new minimum on the ground state surface, leading to a new chemical species 
(product). This may have positive value as a way for synthesizing new mole- 
cules or negative value because of photodegradation of materials or because 
of photochemically-induced cancers. Either way the photochemical reaction 
must occur quickly enough that it can compete with other decay processes. 
Photochemical reactions almost always occur via photochemical funnels where 
excited-state and ground-state surfaces come together, either almost touch- 
ing (avoided crossing) or crossing (conical intersection). These funnels play a 
role in photochemical reactions similar to transition states for thermal reac- 
tions. However it must be kept in mind that these funnels may be far from 
the Frank-Condon region on the excited-state potential energy surface, either 
because there is an easy energetically- "downhill" process or because, unless 
the absorption wavelength can be carefully tuned to a known vertical excita- 
tion energy, the system will typically arrive in an electronically-excited state 
with excess dynamical energy which can be used to move from one excited- 
state potentiall energy surface valley over a transition state to funnels in 
another basin of the excited-state potential energy surface. While conical 
intersections are forbidden in diatomic molecules, they are now believed to 
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be omnipresent in the photochemistry of polyatomic molecules where the 

presence of an avoided crossing in potential energy surface cross-sections ob- 
tained from oversimplified models (e.g., those making simplifying symmetry 
assumptions) usually indicates a nearby conical intersection. A particularly 
striking example is provided by experimental and theoretical evidence that 
the fundamental photochemical reaction involved in vision passes through 
a conical intersection [PoUi 2010]. For these reasons, modern photochemical 
modeling often involves some type of dynamics and, when this is not possible, 
at least focuses on finding conical intersections that can explain the reaction. 

While a single-reference electronic structure method may be adequate for 
describing photophysical processes, the usual standard for describing pho- 
tochemical processes is a multireference electronic structure method such 
as the complete active space self-consistent field (CASSCF) method. (See 
Ref. [Helgaker 2000] for a review of modern quantum chemical methods.) This 
is because the first approximation to the wave function along a reaction path- 
way is as a linear combination of the wave functions of the initial reactants 
and the final products. Since CASSCF is both computationally heavy and 
requires a high-level of user intervention, a simpler method such as TDDFT 
would be very much welcome. Early work in TDDFT in quantum chemistry 
foresaw increasing applications of TDDFT in photochemical modeling. For 
example, avoided crossings between cross-sections of excited-state potential 
energy surfaces may be described with TDDFT because of the multireference- 
likc nature of TDDFT excited states [Casida 1998]. However great attention 
must also be paid to problems arising from the use of approximate function- 
als [Casida 2002]. In particular, the TDDFT Tamm-Dancoff approximation 
(TDA) [Hirata 1999] was found to give improved shapes of excited-state po- 
tential energy smfaccs [Casida 2000, Cordova 2007], albeit at the price of loos- 
ing oscillator strength sum rules. A major advance towards serious investiga- 
tions of TDDFT for describing photoprocesses came with the implementation 
of analytical derivatives for photochemical excited states in many electronic 
structure programs [Van Caillie 1999, Van Caillie 2000, Furche 2002, Hutter 
2003, Rappoport 2005, Doltsinis 2005, Scalmani 2006]. This made it possible 
to relax excited-state geometries and to calculate Stokes shifts within the 
framework of TDDFT. In fact, TDDFT has become a standard part of the 
photochemical modeler's toolbox. It is typically used for calculating absorp- 
tion spectra and exploring excited-state potential energy surfaces around the 
Frank-Condon region. TDDFT also serves as a rapid way to gain the chem- 
ical information needed to carry out subsequent CASSCF calculations. (See 
e.g., Refs. [Diau 2001a, Diau 2001b, Diau 2002, S0lling 2002, Diau 2003] 
for some combined femtosecond spectroscopy/theoretical studies of photo- 
chemical reactions which make good use of TDDFT.) It would be nice to 
bo able to use a single method to model entire photochemical processes. 
The advent of mixed TDDFT/classical surface- hopping Tully-type dynamics 
[Tapavicza 2007, Werner 2008, Tapavicza 2008, Tavernelh 2009, Tavernelh 
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2009, Barbatti 2010] is giving us a way to extend the power of TDDFT to 
the exploration of increasingly complicated photochemical processes. 

The rest of this chapter is organized as follows: The next section reviews 
non-Born-Oppenheimer phenomena from a wave-function point of view, with 
an emphasis on mixed quantum/classical dynamics. This sets the stage for 
our discussion of TDDFT for real-time dynamics and conical intersections in 
Sec. [Ol We sum up in Sec. [TH 

1.2 Wave- Function Theory 

Most likely anyone who has made it this far into this chapter has seen the 
Born-Oppenheimer approximation at least once, if not many times. How- 
ever it is relatively rare to find good discussions that go beyond the Born- 
Oppenheimer approximation [Doltsinis 2002, Cederbaum 2004]. This section 
tries to do just this from a wave-function point of view, in preparation for 
a discussion of TDDFT approaches to the same problems in the following 
section. We first begin by reviewing (again!) the Born-Oppenheimer approxi- 
mation, but this time with the point of view of identifying the missing terms. 
We then discuss mixed quantum/classical approximations, and end with a 
discussion of the pathway method and ways to find and characterize conical 
intersections. We shall use Hartree atomic units (h = = e — 1) though- 
out and adapt the convention in this section that electronic states are labeled 
by small Latin letters, while nuclear degrees of freedom are labeled by capital 
Latin letters. 

Born-Oppenheimer Approximation and Beyond 

As is well-known, the Born-Oppenheimer approximation relies on a sepa- 
ration of time scales: Since electrons are so much lighter and so move so much 
faster than nuclei, the electrons may be thought of as moving in the field of 
nuclei which are "clamped" in place and the nuclei move in a field which 
is determined by the mean field of the electrons. The Born-Oppenheimer 
approximation provides a precise mathematical formulation of this physical 
picture. Our interest here is in where the Born-Oppenheimer approximation 
breaks down and what terms are needed to describe this breakdown. 

Consider a molecule composed of M nuclei and N electrons. Denote 
the nuclear coordinates by R = {Ri, R2, ■ ■ ■ , Rm) and electronic coordi- 
nates by r = (ri, r2, . . . , Tat). The full Hamiltonian, i?(R, r) = T„(R) -|- 
ife(r;R) + Km(R), is the sum of an electronic Hamiltonian, He(r;R) — 
Te{r) + Ven{^', R)+Ke(r), with its electron kinetic energy, Tg, electron-nuclear 
attraction, Ven, and electron-electron repulsion, Vee, with the missing nuclear 
terms — namely the nuclear kinetic energy, r„, and the nuclear- nuclear repul- 
sion, Vnn- Solving the time-dependent Schrodinger equation. 



H{R, r)${R, X, t) = i^<?(R, X, t) , 



(1.1) 
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is a formidable {N + A/)-body problem, (x denotes inclusion of electron spin. 
We have decided to omit nuclear spin, though this should be included in 
principle when identical nuclei with spin are present, such as in the case of 
ortho- and para- hydrogen.) That is why the Born-Oppenheimer expansion 
(which is not yet the Born-Oppenheimer approximation!), 

'?(R,x,0 -^fj(x;R)xj(R,i), (1-2) 

3 

is used, where the electronic wave functions are solutions of the time- independent 
electronic problem in the field of clamped nuclei, i?e(r; R-)^^^^; R) = 
-E|(R)!f'j(x; R). Inserting the Born-Oppenheimer expansion [Eq. ()1.2p ] into 
the full Schrodinger equation [Eq. (jl.ip ]. left multiplying by !f'*(x;R), and 
integrating over x gives the time-dependent Schrodinger equation for the 
nuclear degrees of freedom, 

(t:„(R) + K,(R)) x.(R,0 + E v;j(R)x,(R,i) - »^x.(R,i) • (i-3) 

i 

Here, T^i(R) — EfiTl) + y„„(R), is the adiabatic potential energy surface 
for the ith electronic state. [Notice that this is a different use of the term 
"adiabatic" than in the TDDFT "adiabatic approximation" for the exchange- 
correlation (xc) functional.] The remaining part, Vij (R), is the hopping term 
which couples the Ith and Jth PESs together. As is well known, the Born- 
Oppenheimer approximation neglects the hopping terms, ^r„(R) + Vi(R)^ 

X^{R,t)=^{^/^t)x^{^R,t). 

We, on the other hand, are interested in precisely the terms neglected 
by the Born-Oppenheimer approximation. The hopping term is given by, 

y,,,(R)x,(R,t) = - E/(l/2m/) {g\'](R) + 2F(5(R) • V,) x,(R,t), where, 
G[^j(JV) = J ^*{x-R){Vj<I'j{x;'R)) dx = (i|V^|j), is the scalar coupling 

matrix and, PfJCR) = / <f*(a;;R) (V/*;, (a;; R)) dx ^ (i|V/|j), is the deriva- 
tive coupling matrix [Cederbaum 2004]. Note that the derivative coupling 
matrix is also often denoted ^ and called the nonadiabatic coupling vector 
[Doltsinis 2002] . Here we have introduced a compact notation for some com- 
plicated objects: The scalar coupling matrix is simultaneously a function of 
the nuclear coordinates, a matrix in the electronic degrees of freedom, and 
a vector in the nuclear degrees of freedom, and a matrix in the electronic 
degrees of freedom. The derivative coupling matrix is simultaneously a func- 
tion of the nuclear coordinates, a matrix in the electronic degrees of freedom, 
a vector in the nuclear degrees of freedom and a vector in the three spatial 
coordinates of the Ith nucleus. 

Interestingly the scalar coupling matrix and derivative coupling matrix 
are not independent objects. Rather, making use of the resolution of the 
identity for the electronic states, it is straightforward to show that. 
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2Fi^V(R).Vi. We may then rewrite the time-dependent nuclear equation (|1.3 
as, 



y— 

^ 2mi 



k 




which is known as the group Born-Oppenheimer equation [Cederbaum 2004]. 
Evidently this is an equation which can be solved within a truncated manifold 
of a few electronic statees in order to find fully quantum mechanical solutions 
beyond the Born-Oppenheimer approximation. 

More importantly for present purposes is that Eq. (|1.4p brings out the 
importance of the derivative coupling matrix. The derivative coupling matrix 
can be rewritten as, 

_ (*| (v,i7e(R)) \j) - <5.,, V.Sf (R) 

Since this equation is basically a force-like term, divided by an energy differ- 
ence, we see that we can neglect coupling between adiabatic potential energy 
surfaces when (i) the force on the nuclei is sufhciently small (i.e., the nuclei 
are not moving too quickly) and (ii) when the energy difference between po- 
tential energy surfaces is not too large. These conditions often break down in 
funnel regions of photochemical reactions. There is then a tendency to follow 
diabatic surfaces, which may be defined rigorously by a unitary transforma- 
tion of electronic states (when it exists) to a new representation satisfying the 
condition, i^j-^^(R) w 0. The advantage of the diabatic representation (when 
it exists, which is not always the case) is that it eliminates the off-diagonal 
elements of the derivative coupling matrix in the group Born-Oppenheimer 
equation [Eq. (|1.4p ]. hence eliminating the need to describe surface hopping. 
At a more intuitive level, the character of electronic states tends to be pre- 
served along diabatic surfaces because {i\dj/dt) — R ■ (j|Vj') = R ■ Fij « 
in this representation. For this reason, it is usual to trace diabatic surfaces 
informally in funnel regions by analyzing electronic state character, rather 
than seeking to minimize the nonadiabatic coupling vector. Avoided cross- 
ings of adiabatic surfaces are then described as due to configuration mixing 
of electronic configurations belonging to different diabatic surfaces. 

Mixed Quantum/ Classical Dynamics 

Solving the fully quantum-mechanical dynamics problem of coupled elec- 
trons and nuclei is a challenge for small molecules and intractable for larger 
molecules. Instead it is usual to use mixed quantum/classical methods in 
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which the nuclei are described by Newtonian classical mechanics while the 
electrons are described by quantum mechanics. Dividing any quantum sys- 
tem into two parts and then approximating one using classical mechanics 
is the subject of on-going research [Kapral 2006]. In general, no rigorous 
derivation is possible and wave- function phase information (e.g., the Berry 
phase) is lost which may be important in some instances. Nevertheless mixed 
quantum/classical approximations are intuitive: Most nuclei (except perhaps 
hydrogen) are heavy enough that tunneling and other quantum mechanical 
effects are minor, so that classical dynamics is often an a priori reasonable 
first approximation. Of course, rather than thinking of a single classical tra- 
jectory for the nuclear degree of freedom, we must expect to think in terms of 
ensembles (or "swarms" ) of trajectories which are built to incorporate either 
finite temperature effects or to try to represent quantum mechanical prob- 
ability distributions or both. The purpose of this subsection is to introduce 
some common mixed quantum/classical methods. 

The most elementary mixed quantum/classical approximation is Ehren- 
fest dynamics. According to Ehrenfest's theorem [Ehrenfest 1927], Newton's 
equations are satisfied for mean values in quantum systems, d{r)/dt = {p)/m 
and d{p)/dt = —{VV) . Identifying the position of the nuclei with their mean 
value, we can then write an equation, mjRi{t) = — V/t^(R(t)), whose phys- 
ical interpretation is that the nuclei are moving in the mean field of the 
electrons. Here 

y(R(t)) = {<FiK,t)\H,{Ii{tm{Ii,t)) +Vnnmt)) , (1-6) 

where the electronic wave function is found by solving the time-dependent 
equation, -ffe(x, R(t))^'(x; R, t) = i{d/dt)\P{-K; R, t). While Ehrenfest dynam- 
ics has been widely and often successfully applied, it suffers from some im- 
portant drawbacks. The first drawback is that the nuclei always move on 
average potential energy surfaces, rather than adiabatic or diabatic surfaces, 
even when far from fimnel regions where the nuclei would be expected to 
move on the surface of a single electronic state. While this is serious enough, 
since it suggests errors in calculating branching ratios (i.e., relative yields of 
different products in a photorcaction), a more serious drawback is a loss of 
microscopic reversibility. That is, the temporal variation of the mean poten- 
tial energy surface depends upon past history and can easily be different for 
forward and reverse processes. 

A very much improved scheme is the fewest switches method of Tully 
[Tully 1990, Hammes-Schiffer 1994]. Here the nuclei move on well-defined adi- 
abatic potential energy surfaces, miRi{t) = — V/Vi(R(t)), and the electrons 
move in the field of the moving nuclei, ffe(r; R(t))!f'(r, t) = i{d/ dt)^{v,t). 
To determine the probability that a classical trajectory describing nuclear 
motion hops from one electronic potential energy surface to another, we ex- 
pand SI/(T,t) = ^yjj !f'm(r; R(t))Cm(t), in solutions of the time-independent 
Schrodinger equation, .H"(r; R(t))^'„(r; R(t)) = £'„(R(t))!2>„(r; R(t)). The 
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probability of finding the system on surface m is then given by, Pm{t) = 

ICmCi)!^- The coefficients may be obtained in a dynamics calculation by in- 
tegrating the first-order equation, Cm{t) = —iEm{t)Cm{t) — 
X Cn{t). A not unimportant detail is that the nonadiabatic coupling elements 
need not be calculated explicitly, but instead can be calculated using the fi- 
nite difference formula, (m{t + At/2)\h{t + At/2)) = 

[{m{t)\n{t + St)) — {m{t + At)\n{t))] /{2At). In practice, it is also important 
to minimize the number of surface hops or switches in order to keep the cost 
of the dynamics calculation manageable. Tully accomplished this by introduc- 
ing his fewest-switches algorithm which is a type of Monte Carlo procedure 
designed to correctly populate the different PESs with a minimum of surface 
hopping. Briefly, the probability of jumping from surface m to surface n in 
the interval {t, t + At) is given by gm->-n{t, At) = Pm,n{t)At/Pm,m{t) where 
Pm,n{t) = Cm{t)C*{t). A random number ^ is generated with uniform prob- 
ability on the interval (0, 1) and compared with gm^n{t, At). The transition 

m ^ n occurs only if Pi™-'^ < ,e < Pn^^ where Pi™^ = E;=i ™ Pn,i is the 
sum of the transition probabilities for the first m states. Additional details 
of the algorithm, beyond the scope of this chapter, involve readjustment of 
nuclear kinetic energies and the fineness of the numerical integration grid for 
the electronic part of the calculation with respect to that of the grid for the 
nuclear degrees of freedom. 

It is occasionally useful to have a simpler theory for calculating the prob- 
ability of potential energy surface hops which depends only on the potential 
energy surfaces and not on the wave functions. Such a theory was suggested by 
Landau [Landau 1932] and Zener [Zener 1932] (see also Wittig [Wittig 2005]). 
Their work predates the modern appreciation of the importance of conical 
intersections and so focused on surface hopping at avoided crossings. The 
Landau-Zener model assumes that surface hopping occurs only on the surface 
where the two diabatic surfaces cross that give rise to the avoided crossing 
where the surface hopping occurs. After some linearizations and an asymp- 
totic limit, it is possible to arrive at a very simple final formula, 

^ = ^''^{~h{d\AEaS/dt)) ' ^^-^^ 

for the probability of hopping between two potential energy surfaces. This 
formula is to be applied at the point of closest approach of the two potential 
energy surfaces where the energy difference is AEadia- However d\AEciia\/dt 
is evaluated as the maximum of the rate of change of the adiabatic energy 
difference as the avoided crossing is approached. While not intended to be 
applied to conical intersections, it is still quite applicable in photodynamics 
calculations since trajectories rarely go exactly through a conical intersection. 

Pathway Method 

Dynamics calculations provide a swarm of reaction trajectories. The "path- 
way method" provides an alternative when dynamics calculations are too ex- 
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pensive or a simplified picture is otherwise desired, say, for interpretational 
reasons. The pathway method consists of mapping out minimum energy path- 
ways between the initial Franck-Condon points obtained by vertical excita- 
tions and excited-state minima or conical intersections. Although analogous 
to the usual way of finding thermal reaction paths, it is less likely to be a 
realistic representation of true photoprocesses except in the limit of threshold 
excitation energies since excess energy is often enough to open up alternative 
pathways over excited-state transition states. While the necessary ingredi- 
ents for the photochemical pathway method are similar to those for thermal 
reactions, conical intersections are a new feature which is quite different from 
a thermal transition state. This section provides a brief review for finding 
and characterizing conical intersections. 

The notion of a conical intersection arises from a relatively simple ar- 
gument [Yarkony 2001]. The potential energy surface of a molecule with / 
internal degrees of freedom is an /-dimensional hypersurface in an (/ -I- 1)- 
dimensional space (the extra dimension is the energy axis). If two potential 
energy surfaces simply cross "without seeing each other" , then the crossing 
space is characterized by the constraint. 



makes the crossing space (/— l)-dimensional. However in quantum mechanics, 
we also have the additional constraint. 



This makes the crossing space (/ — 2)-dimcnsional. This means that there 
will be two independent directions in hyperspace in which the two potential 
energy surfaces will separate. These two directions define a branching plane. 
Within the 3-dimensional space defined by the energy and the branching 
plane, the conical intersection appears to be a double cone (Fig. II. 6p . the 
point of which represents an entire (/ — l)-dimensional space. Of course, 
/ = 1 for a diatomic and no conical intersection is possible. This is the origin 
of the well-known avoided crossing rule for diatomics. Here we are intererested 
in larger molecules where the low dimensionality of the branching space in 
comparison with the dimensionality of the parent hyperspace can make the 
conical intersection hard to locate and characterize. 

In the pathway method, the system simply goes energetically downhill un- 
til two potential energy surfaces have the same energy [Eq. (|1.8|) ]. The resul- 
tant intersection space must be analyzed and the branching plane extracted 
so that the surface crossing region can be properly visualized and interpret- 
ted. In order to do so, let us recall a result from first-year calculus. Imagine 
a trajectory, R(t), depending upon some parameter r within the conical 
intersection surface. Then VC(R) must be perpendicular to the conical in- 
tersection for any constraint function C(R) — because, = (iC(R(r))/dr = 



(1.8) 



(1.9) 
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VC(R) ■ (dR/dr). and we can always choose dR/dr ^ 0. Taking the gra- 
dient of Eq. (|1.9p defines the derivative coupling vector, ^ — Vi/i.j(R), 
while taking the gradient of Eq. (|1.8p defines the gradient difference vector, 
= Vi?i(R) — Vi?j(R). Together the derivative coupling vector and gra- 
dient difference vector are referred to as the branching vectors which char- 
acterize the branching plane. [Note that the derivative coupling vector is 
essentially the numerator of the derivative coupling matrix expression given 
in Eq. (|1.5|) . This confusion of nomenclature is unfortunate but present in 
the literature.] 

These branching plane conditions are needed as constraints in the ex- 
ploration of the conical intersection hyperspace when seeking the minimum 
energy conical intersection or the first-order saddle point in conical inter- 
seection. In particular, there has been considerable effort devoted to the 
problem of developing efhcient algorithms for finding minimum energy points 
within the conical intersection space [Koga 1985, Atchity 1991, Ragazos 1992, 
Yarkony 1996, Domcke 1997, Izzo 2000]. Furthermore, an automated system- 
atic exploration method for finding minimum energy conical intersections has 
very recently developed [Maeda 2009] . First-order saddle points and the cor- 
responding minimum energy pathways both within the conical intersection 
hypersurface have been proposed to be important in dynamical trajectory 
simulations, and an optimization method was developed for such high-energy 
points within the conical intersection hypersurface [Sicilia 2008] . Some of the 
minimum energy conical intersection optimizers use the branching plane con- 
ditions explicitly to keep the degeneracy of the two adiabatic states during 
optimizations [Manaa 1993, Bearpark 1994, Anglada 1997], making explicit 
use of both the derivative coupling vector and gradient difference vector at ev- 
ery step. Most well-estabilished optimization algorithms assume smoothness 
of the function to be optimized. Since the potential energy surface necessarily 
has a discontinuous first derivative in the vicinity of a conical intersection, 
the above-mentioned algorithms for finding minimum energy conical inter- 
sections have required access to the gradient difference vector and deriva- 
tive coupling vectors. The gradient difference vector can easily be obtained 
from analytical gradients, if available, or by numerical energy differentiation 
if analytical gradients are not yet available. However methods for finding 
the derivative coupling vector are not yet available for all methods since 
implementation of an analytical derivative method is often regarded as a 
prerequisite [Ciminelli 2004, Maeda 2010]. Some approaches make use of a 
penalty function to get around the need to calculate the derivative coupling 
vector and these have proven very useful for finding minimum energy coni- 
cal intersection regions without the need for the derivative coupling vector 
[Levine 2008] . This is especially important for methods such as renormalized 
coupled-cluster theories and TDDFT or free-energy methods for which the 
electronic wave function is not completely defined, considerably complicating 
the problem of how to calculate derivative coupling vector matrix elements. 
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However, convergence of penalty function methods is in general slower than 
methods which make explicit use of the branching plane constraints, espe- 
cially if tight optimization of the energy difference, {Ei — Ej), is desired 
[Keal 2007]. 

1.3 TDDFT 

The last section discussed the basic theory of non-Born-Oppenheimer dy- 
namics and conical intersections from a wave-function point of view. We now 
wish to see to what extent we can replace wave-function theory with what 
we hope will be a simpler DFT approach. As usual in DFT, we seek both the 
guiding light of formal rigor and pragmatic approximations that work. We 
will take a more or less historical approach to presenting this material. In 
this section, upper case Latin indices designate electronic states, while lower 
case Latin indices designate orbitals. 

One of the early objectives of TDDFT was to allow simulations of the 
behavior of atoms and clusters in intense laser fields, well beyond the linear- 
response regime and too complex to be handled by comparable wave- function 
methods. The closely related topic of ion-cluster collisions was studied early 
on using TDDFT in a very simplified form [Yabana 1998]. The Ehrenfest 
method was the method of choice for TDDFT simulations coupling elec- 
tronic and nuclear degrees of freedom in this area. The main problem is how 
to take the expectation value in Eq. (jl.6p . This is solved pragmatically by 
using, V[p{t)]{'R.{t)) =T,[^,{t)] + J {v,n{r,-R{t)) + vappiir,t)) piv,t) dr + 
(1/2) / / {p{ri,t)p{r2,t)/ri2)dridr2+E,,[p{t)] + Vnn{Ii{t)), where T, is the 
usual Kohn-Sham noninteracting kinetic energy evaluated using the (now 
time-dependent) Kohn-Sham determinant Ven is the electron- nuclear at- 
traction potential, Vappi is the potential for any applied electric field, and 
Exc is the usual Hohenberg-Kohn-Sham xc-energy. Note that the presence 
of Exc is highly reminiscent of the TDDFT adiabatic approximation that 
the xc-potential should react instantaneously and without memory to any 
temporal change in the time-dependent charge density. Among the notable 
work done with this approximation is early studies of the dynamics of sodium 
clusters in intense laser fields [Calvayrac 1998], the development of the time- 
dependent electron localization function [Burnus 2005], and (more recently) 
the study of electron-ion dynamics in molecules under intense laser pulses 
[Kawashita 2009] . Besides limitations associated with the TDDFT adiabatic 
approximation, the TDDFT Ehrenfest method suffers from the same intrin- 
sic problems as its wave-function brother — namely that it is implicitly based 
on an average potential energy surface and so does not provide state-specific 
information, and also suffers from problems with microscopic irreversibility. 

To our knowledge, the first DFT dynamics on a well-defined excited-state 
potential energy surface was not based upon TDDFT but rather on the older 
multiplet sum method of Ziegler, Rank, and Baerends [Ziegler 1977, Daul 
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1994]. This was the work of restricted open-shell Kohn-Sham (ROKS) for- 
malism of Irmgard Frank et al. [Frank 1997] who carried out Carr-Parinello 
dynamics for the open-shell singlet excited state ^(j,a) using the multiplet 
sum method energy expression, Es = 2E[$'^J] - where 'Pf; is the 

Kohn-Sham determinant with the ia spin-orbital replaced with the ar spin- 
orbital. Such a formalism suffers from all the formal difficulties of the multi- 
plet sum method, namely that it is just a first-order estimate of the energy 
using a symmetry-motivated zero-order guess for the excited-state wave func- 
tion and assumes that DFT works best for states which are well-described 
by single determinants. Nevertheless appropriate use of the multiplet sum 
method can yield results similar to TDDFT. A recent application of this 
method is to the study of the mechanism of the electrocylic ring opening of 
diphenyloxirane [Friedrichs 2009]. 

The implementation of TDDFT excited-state derivatives in a wide variety 
of programs not only means that excited-state geometry optimizations may 
be implemented, allowing the calculation of the Stokes shift between absorp- 
tion and fluorescence spectra, but that the pathway method can be imple- 
mented to search for conical intersections in TDDFT. Unless nonadiabatic 
coupling matrix elements can be calculated within TDDFT [vide infra), then 
a penalty method should be employed as described in the previous section 
under the pathway method. This has been done by Levine, Ko, Quenneville, 
and Martinez using conventional TDDFT [Levine 2006] and by Minezawa 
and Gordon using spin- flip TDDFT [Minezawa 2009]. We will come back to 
these calculations later in this section. 

Thee most recent approach to DFT dynamics on a well-defined excited- 
state potential energy surface is Tully-type dynamics [Tully 1990, Hammes- 
Schiffer 1994, Tully 1998] applied within a mixed TDDFT/classical trajectory 
surface-hopping approach. Surface- 
hopping probabilities can be calculated from potential energy surfaces alone 
within the Landau-Zener method [Eq. ()1.7p ]. however a strict application 
of Tully 's method requires nonadiabatic coupling matrix elements as input. 
Thus a key problem to be addressed is how to calculate nonadiabatic coupling 
matrix elements within TDDFT. Initial work by Craig, Duncan, and Prezhdo 
used a simple approximation which neglected the xc-kernel [Craig 2005]. A 
further approximation, criticized by Maitra [Maitra 2006], has been made 
by Craig and co-workers [Craig 2005, Habenicht 2006] who treated the elec- 
tronic states as determinants of Kohn-Sham orbitals which are propagated 
according to the time-dependennt Kohn-Sham equation. This means that 
neither the excitation energies nor the associated forces could be consid- 
ered to be accurate. The first complete mixed TDDFT/classical trajectory 
surface-hopping photodynamics method was proposed and implemented by 
Tapavicza, Tavernelli, and Rothlisberger [Tapavicza 2007] in a development 
version of the CPMD code. It was proposed that the nonadiabatic coupling 
matrix elements be evaluated within Casida's ansatz [Casida 1995] which 
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was originally intended to aid with the problem of assigning excited states 
by considering a specific functional form for an approximate excited-state 
wave function. For the TDA, the Casida ansatz takes the familiar form, 
= '^iaa'^'ia ■^io.(T- In fact, matrix elements between ground and excited 
states may be calculated exactly in a Casida-like formalism because of the 
response theory nature of Eq. dUS]) [Chernyak 2000, Hu 2007, Send 2010]. 
Test results show reasonable accuracy for nonadiabatic coupling matrix ele- 
ments as long as conical intersections are not approached too closely [Baer 
2002, Hu 2007, Tavernelh 2009, Tavernelli 2009, Send 2010]. One likely rea- 
son for this is the divergence of Eq. (|1.5|) when Ei = Ej. Hu and Sugino 
attempted to further improve the accuracy of nonadiabatic coupling matrix 
elements by using average excitation energies [Hu 2007]. The problem of cal- 
culating nonadiabatic coupling matrix elements between two excited states 
is an open problem in TDDFT, though the ability to calculated excited- 
state densities [Furche 2002] suggests that such matrix elements could be 
calculated from double response theory using Eq. ()1.5p . Soon after the im- 
plementation of mixed TDDFT/classical trajectory surface- hopping photo- 
dynamics in CPMD, a very similar method was implemented in TurboMol 
and applied [Werner 2008, Miric 2008, Barbatti 2010]. A version of Turbo- 
Mol capable of doing mixed TDDFT/classical trajectory surface- hopping 
photodynamics using analytic nonadiabatic coupling matrix elements has re- 
cently appeared [Send 2010] and has been applied used to study the photo- 
chemistry of vitamin-D [Tapavicza 2011]. Time-dependent density- functional 
tight-binding (TDDFTB) may be regarded as the next step in a multiscale 
approach to the photodynamics of larger systems. From this point of view, 
it is interesting to note that mixed TDDFTB/classical trajectory surface- 
hopping photodynamics is also a reality [Mitric 2009] . Given the increasingly 
wide-spread nature of implementations of mixed TDDFT classical trajectory 
surface-hopping photodynamics, we can only expect the method to be in- 
creasingly available to and used by the global community of computational 
chemists. 

Before going further, let us illustrate the state-of-the-art for TDDFT when 
applied to non-Born-Oppenheimer dynamics and conical intersections. We 
will take the example of the photochemical ring opening of oxirane (struc- 
ture I in Fig. II. 2|) . While this is not the "sexy application" modeling of some 
biochemical photoprocess, the photochemistry of oxiranes is not unimportant 
in synthetic photochemistry and, above all, this is a molecule where it was 
feh that TDDFT "ought to work" [Cordova 2007]. A first study showed that 
a main obstacle to photodynamics is the presence of triplet and near singlet 
instabilities which lead to highly underestimated and even imaginary excita- 
tion energies as funnel regions are approached. This is illustrated in Fig. 11.31 
for C2v ring opening. While the real photochemical process involves asymet- 
ric CO ring-opening rather than the symmetric C2v CC ring-opening, results 
for the symmetric pathway have the advantage of being easier to analyze. 



1 Real-Time Dynamics and Conical Intersections 



15 



/\ 

H-C— C-H 
I I 
H H 



(1) 



:0: 
I . 
H-C— C- 
I I 
H H 





- I © 
H-C— C-H 
I I 
H H 



(2) 



H-C- 



(4) 



H 

■C-H 
H 



0- 

II 

H-C- 



H 

■C-H 
H 



(3) 



Fig. 1.2. Mechanism proposed by Gomer and Noyes in 1950 for the photochemical 
ring opening of oxirane. Reprinted with permission from E. Tapvicza, I. Tavernelli, 
U. Rothlisberger, C. Filippi, and M. E. Casida, J. Chem. Phys. 129, 124108 (2008). 
Copyright 2008, American Institute of Physics. 



The figure shows that applying the TDA strongly attenuates the instability 
problem, putting most curves in the right energy range. Perhaps the best way 
to understand this is to realize that, whereas time-dependent Hartree-Fock 
(TDHF), is a nonvariational method and hence allows variational collapse of 
excited states, TDA TDHF is the same as configuration interaction singles 
(CIS) which is variational. There is however still a cusp in the ground state 
curve as the ground state configuration changes from to (cr*)^. Accord- 
ing to a traditional wave-function picture, these two states, which are each 
double excitations relative to each other should be included in configuration 
mixing in order to obtain a proper description of the ground state potential 
energy surface in the funnel region [Cordova 2007, Huix-Rotllant 2010]. 

Figure [L4l shows an example of mixed TDA TDPBE/classical trajectory 
surface-hopping calculations for the photochemical ring-opening of oxirane 
with the initial photoexcitation prepared in the ^(n,3pz) state. Part (b) of 
the figure clearly shows that more than one potential energy surface is popu- 
lated after about 10 fs. The Landau-Zener process is typical of the dominant 
physical process which involves an excitation from the HOMO nonbonding 
lone pair on the oxygen initially to a 3pz Rydberg orbital. As the reaction 
proceeds, the ring opens and the target Rydberg orbital rapidly changes char- 
acter to become a CO a* antibonding orbital (Fig. II. 5| ). Actual calculations 
were run on a swarm of 30 trajectories, confirming the mechanism previously 
proposed Gomer-Noyes mechanism [Gomer 1950] (Fig. ll.2( ). but also confirm- 
ing other experimental by-products and giving unprecedented state-specific 
reaction details such as the orbital description briefly described above. 

The oxirane photochemical ring-opening passes through a conical inter- 
section, providing a concrete example of a conical intersection to study with 
TDDFT. We now return to the study by Levine, Ko, Quenneville, and Mar- 
tinez of conical intersections using conventional TDDFT [Levine 2006] who 
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Fig. 1.3. Comparison of TDA TDLDA and diffusion Monte Carlo curves for C2v 
ring opening of oxirane. Reprinted with permission from F. Cordova, L. Joubert 
Doriol, A. Ipatov, M. E. Casida, C. Filippi, and A. Vela, J. Chem. Phys. 127, 
164111 (2007). Copyright 2007, American Institute of Physics. 



noted that strict conical intersections are forbidden by the TDDFT adiabatic 
approximation for the simple reason that there is no coupling matrix ele- 
ment [Eg. II. 9j to zero out between the ground and excited states. Figure [T^ 
shows a CASSCF conical intersection close to the oxirane photochemical 
funnel. Also shown are the TDA TDDFT surfaces calculated with the same 
CASSCF branching coordinates. Interestingly the CASSCF and TDDFT con- 
ical intersections look remarkably similar. However closer examination shows 
that the TDDFT "conical intersection" is actually two intersecting cones 
rather than a true conical intersection, confirming the observation of Levine 
et al. This was analyzed in detail in Ref. [Tapavicza 2008] where it was con- 
cluded that the problem is that we are encountering effective noninteracting 
f-representability. True noninteracting u-representability means that there is 
no noninteracting system whose ground state gives the ground state density 
of the interacting system. This only means that there is some excited state 
of the noninteracting system with integer occupation number which gives the 
ground state density of the interacting system. What we call effective nonin- 
teracting u-representability is when the LUMO falls below the HOMO (or, in 
the language of solid-state physics, there is a "hole below the Fermi level"). 
This is exactly what frequently happens in the funnel region. 

Spin-flip (SF) TDDFT [Slipchenko 2003, Shao 2003, Wang 2004] offers 
one way to circumvent some of the problems of effective noninteracting v- 
representability in funnel regions. This is because we can start from the low- 
est triplet state which has fewer effective noninteracting u-representability 
problems and then use SFs to obtain both the ground state and a doubly- 
excited state. Analytic derivatives are now available for some types of SF- 



1 Real-Time Dynamics and Conical Intersections 



17 




Fig. 1.4. (a) Cut of potential energy surfaces along reaction path of an Landau- 
Zener (dashed line) and a fewest-switches (solid line) trajectory (black, 5*0; blue, 
5*1; green, 5*2; magenta, S3). Both trajectories were started by excitation into the 
^(n, 3pz) state, with the same geometry and same initial nuclear velocities. The 
running states of the Landau-Zener and the fewest-switches trajectory are indicated 
by the red crosses and circles, respectively. The geometries of the Landau-Zener 
trajectory are shown at time a) 0, b) 10, and c) 30 fs. (b) State populations (black, 
So; blue, ^i; green, S2; magenta, 5*3) as a function of the fewest-switches trajectory 
in (a). Reprinted with permission from E. Tapvicza, 1. Tavernelli, U. Rothlisberger, 
C. Filippi, and M. E. Casida, J. Chem. Phys. 129, 124108 (2008). Copyright 2008, 
American Institute of Physics. 



TDDFT [Scth 2010]. Figure O shows that SF-TDDFT works fairly well for 
treating the avoided crossing in the C21, ring-opening pathway of oxirane. 
Minezawa and Gordon also used SF-TDDFT to locate a conical intersec- 
tion in ethylene [Minezawa 2009]. However Huix-Rotllant, Natrajan, Ipatov, 
Wawire, Deutsch, and Casida found that, although SF-TDDFT does give a 
true conical intersection in the photochemical ring opening of oxirane, the 
funnel is significantly shifted from the position of the CASSCF conical inter- 
section [Huix-Rotllant 2010]. The reason is that the key funnel region involves 
an active space of over two orbitals which is too large to be described accu- 
rately by SF-TDDFT. 
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Fig. 1.5. Change of character of the active state along the reactive Landau-Zener 
trajectory, shown in Fig. 11.41 Snapshots were taken at times (a) 2.6, (b) 7.4, (c) 12.2, 
and (d) 19.4 fs. For (a) and (b), the running state is characerized by a transition from 
the highest occupied molecular orbital (HOMO) to the lowest unoccupied molecular 
orbital (LUMO) plus one (LUMO+1), while for (c) and (d) it is characterized by 
a HOMO-LUMO transition due to orbital crossing. Note that the HOMO remains 
the same oxygen nonbonding orbital throughout the simulation. Reprinted with 
permission from E. Tapvicza, 1. Tavernelli, U. Rothlisberger, C. Filippi, and M. E. 
Casida, J. Chem. Phys. 129, 124108 (2008). Copyright 2008, American Institute of 
Physics. 
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Fig. 1.6. Comparison of the So and Si potential energy surfaces calculated using 
different methods for the CASSCF branching coordinate space. M. Huix-RotUant, 
B. Natarajan, A. Ipatov, C. M. Wawire, T. Deutsch, and M. E. Casida, Phys. Chem. 
Chem. Phys. 12, 12811 (2010) — Reproduced by permission of the PCCP Owner 
Societies. 
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Fig. 1.7. C2v potential energy curves: full calculation (solid lines), two-orbital 
model (dashed lines). M. Huix-RotUant, B. Natarajan, A. Ipatov, C. M. Wawire, 
T. Deutsch, and M. E. Casida, Phys. Chem. Chem. Phys. 12, 12811 (2010) — 
Reproduced by permission of the PCCP Owner Societies. 

There are other ways to try to build two- and higher-excitation char- 
acter into a DFT treatment of excited states. Let us mention here only 
multireference configuration interaction (MRCI)/DFT [Grimme 1999], con- 
strained density functional theory-configuration interaction (CDFT-CI) [Wu 
2007], and mixed TDDFT/many-body theory methods based upon the Bethe- 
Salpeter equation [Romaniello 2009] or the related polarization propagator 
approach [Casida 2005, Huix-RotUant 2010] or the simpler dressed TDDFT 
[Maitra 2004, Cave 2004, Gritsenko 2009, Mazur 2009, Mazur 2010, Huix- 
RotUant 2011] approach. 

All of these may have the potential to improve the DFT-based description 
of funnel regions in photochemical reactions. Here however we must be aware 
that we may be in the process of building a theory which is less automatic and 
requires the high amount of user intervention typical of present day CASSCF 
calculations. This is certainly the case with CDFT-CI which has already 
achieved some success in describing conical intersections [Kaduk 2010]. 

1.4 Perspectives 

Perhaps the essence of dynamics can be captured in a simple sentence: "You 
should from whence you are coming and to where you are going." Of course 
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this rather deterministic statement must be interpreted differently in classi- 
cal and quantum mechanics. Here however we would like to think about its 
meaning in terms of the development of DFT for applications in photopro- 
cesses. Theoretical developments in this area have been remarkable in recent 
years, opening up the possibility for a more detailed understanding of fem- 
tosecond (and now also attosecond) spectroscopy. In this chapter we have 
tried to discuss the past, the present, and a bit of the future. 

The past treated here has been the vast area of static investigation and 
dynamic simulations of photophysical and photochemical processes. We have 
first described more traditonal wave-function techniques. We have also men- 
tionned and made appropriate references to important work on early DFT 
work involving Ehrenfest TDDFT and restricted open-shell Kohn-Sham DFT 
dynamics. Our emphasis has been on photochemical processes involving sev- 
eral PESs, partly because of our own personal experiences, but also because 
photo c/iemica^ processes start out as photophysical processes in the Franck- 
Condon region and then rapidly become more complicated to handle. 

The present-day status of DFT photodynamics is perhaps best repre- 
sented by the recent availability of mixed TDDFT and TDDFTB/classical 
surface-hopping dynamics codes as well as serious efforts to investigate and 
improve the quality of the TDDFT description of photochemical funnel re- 
gions. First applications have already shown the utility of this theory and we 
feel sure that other applications will follow as programs are made broadly 
available to computational scientists. Finally we have ended the last section 
with some speculations about the future concerning the need for explicit 
double- and higher-excitations to correctly describe funnel regions. 

As expected, we could not treat everything of relevance to the chap- 
ter title. Roi Baer's recent work indicating that Berry phase information 
is somehow included in the ground-state charge density is most intriguing 
[Baer 2010]. Also on-going work on multicomponent DFT capable of treat- 
ing electrons and nuclei on more or less the same footing [Kreibich 2001, 
Kreibich 2008] would seem to open up new possibililtics for developing use- 
ful non-Born-Oppcnhcimcr approximations within a DFT framework. We are 
sure that still other potentially relevant work has been unfortunately omitted 
either because of space limitations or for other reasons. 

Do we know where this field is going? Certainly non-Born-Oppenheimer 
photodynamics using some form of DFT is currently a hot and rapidly evolv- 
ing area. Exactly what lies in store may not yet be clear, but what we do 
know is that we are going to have fun getting there! 
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